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We present a Wannier-based method to calculate the Chern-Simons orbital magnetoelectric cou- 
pling in the framework of first-principles density-functional theory. In view of recent developments 
in connection with strong Z2 topological insulators, we anticipate that the Chern-Simons contri- 
bution to the magnetoelectric coupling could, in special cases, be as large or larger than the total 
magnetoelectric coupling in known magnetoelectrics like Cr2 03. The results of our calculations for 
the ordinary magnetoelectrics Cr2 03, BiFe03 and GdA103 confirm that the Chern-Simons contri- 
bution is quite small in these cases. On the other hand, we show that if the spatial inversion and 
time-reversal symmetries of the Z2 topological insulator Bi2Se3 are broken by hand, large induced 
changes appear in the Chern-Simons magnetoelectric coupling. 
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I. INTRODUCTION 

In recent years there has been a significant revival of 
interest in magnetoelectric effects in solids, as surveyed in 
several reviews. 1-4 Potential applications of these materi- 
als have long been discussed 5 ' 6 in areas ranging from the 
optical manipulation and frequency conversion to magne- 
toelectric memories. Of the various quantities that can be 
discussed, the linear magnetoelectric coupling tensor ay 
is clearly of primary interest, as it quantifies the leading- 
order term in the coupling at small fields. We define it 
as 

where Vi is the electric polarization induced by the mag- 
netic field Bj, or equivalcntly, Mj is the magnetization 
induced by the electric field £j. We use SI units (see 
Sec. II A) and the derivatives are to be evaluated at zero 
electric and magnetic field. In the special case that the 
induced response ("P or M) remains parallel to the ap- 
plied field (B or £ ), the tensor a. is purely diagonal with 
equal diagonal elements, and its strength can be mea- 
sured by a dimensionless scalar parameter defined via 

Be 2 

< = 2^h 6ij - (2) 

More generally, depending on the magnetic point group 
of the crystal, ay can have distinct diagonal components 
as well as non-zero off-diagonal ones. 

The linear magnetoelectric response ay can be decom- 
posed into two contributions coming from purely elec- 
tronic and from ionic responses respectively The former 
is defined as the magnetoelectric response that occurs 
when atoms are not allowed to displace in response to the 
applied field, while the latter is defined as the remaining 
lattice-mediated response. One generally expects ionic 
effects to dominate over electronic responses, as for ex- 
ample was shown recently in Ref. 7 and 8 for the case 



of Cr203. Moreover, each of these components can be 
decomposed further into spin and orbital parts, since the 
magnetization induced by the electric field can be decom- 
posed in that way. Here one would naively expect that 
the spin contribution will dominate with respect to the 
orbital one, since orbital moments are usually strongly 
quenched by crystal fields. Mostly for this reason, real- 
istic theoretical calculations of magnetoelectric coupling 
have been developed only for the spin component. 

As shown in Rcfs. 10 and 11 using two complemen- 
tary approaches, the orbital magnetoelectric polarizabil- 
ity (OMP), defined as the contribution of orbital currents 
to the magnetoelectric coupling ay, can be written as 
the sum of three gauge-invariant contributions. One of 
these, first discussed by Qi et al. 12 and Essin et al., 13 
is the Chern-Simons term (CSOMP). Since this contri- 
bution is purely isotropic, we can measure its strength 
using the single parameter 8 as in Eq. (2). In this paper 
we will focus mostly on the CSOMP component of ay. 
From an implementation viewpoint, the CSOMP compo- 
nent is quite different from the other two components of 
the OMP: it can be calculated from a knowledge of the 
ground-state electron wavefunctions alone, but only after 
careful attention is given to the need to choose a smooth 
gauge in discrctized fc-space. 

One of the motivations for the current work is the pos- 
sibility of finding a material whose CSOMP component of 
the linear magnetoelectric tensor will be large compared 
to the total coupling in known magnetoelectric materials. 
As elaborated in more detail in Sec. II, the basis for this 
possibility arises from the before-mentioned theoretical 
developments 14 and the experimental verification of the 
existence of Z2 topological insulators such as Bii-^Sba;, 
Bi2Se3, Bi2Te3 and Sb2Te3. 15-17 Roughly speaking, we 
seek a material that is similar to a Z2 topological insula- 
tor, but having broken inversion and time-reversal sym- 
metries. In order to take the first steps toward search- 
ing for such materials, we have set out to calculate the 
CSOMP component of the magnetoelectric tensor in sev- 
eral compounds of interest using density-functional the- 
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ory. 

The paper is organized as follows. In Sec. II we pro- 
vide theoretical background by reviewing the previously- 
derived 10,11 expression for the a. tensor, and by dis- 
cussing the connection between bulk and surface prop- 
erties in a way that is analogous to the theory of surface 
charge and bulk electric polarization. We also review the 
connection to Z2 topological insulators and make some 
general comments about symmetry. In Sec. Ill we dis- 
cuss the gauge-fixing issues that arise when discretizing 
the CSOMP expression on a fc-point mesh, and show how 
these can be resolved using Wannier-based methods. By 
this route, we arrive at an explicit expression for the 
CSOMP in terms of position matrix elements between 
Wannier functions. We evaluate this expression in the 
density-functional context for several materials of inter- 
est in Sec. IV. Finally, we summarize and give an outlook 
in Sec. V. 



II. BACKGROUND AND MOTIVATION 

In this section we briefly summarize previous work 
from Rcfs. 10 and 11 on the orbital magnetoelectric 
coupling (OMP), describe relationships between bulk 
and surface properties, discuss motivations for this work 
based on the discovery of strong Z2 topological insula- 
tors, and present a brief symmetry analysis. 



A. Units and conventions 

In this paper we use SI units and define a according 
to Eq. (1) using independent field variables £ and B. It 
follows that a has the same units as the vacuum admit- 
tance 1/c/Jo- 18 While this is convenient from the point of 
view of first-principles theory, where B is fixed to zero 
in practice, the more conventional definition in the liter- 
ature is in terms of fixed £ and H fields, in which case 
one has 
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(3) 



and a has units of inverse velocity. In the typical 
case that the magnetic susceptibility of the material is 
negligible, these are related by a EH = afi , and one can 
define a reduced (dimensionless) quantity a r = c^a = 
ca EH . 18 Defined in this way, a r is numerically equal to 
the value of the magnetoelectric coupling in Gaussian 
units using the conventions of Rivera, 19 which in turn 
corresponds to the notation "g.u." ("Gaussian units") in 
some recent papers. 7 ' 9 Furthermore, using the notation 
of Eq. (2) for the isotropic magnetoelectric coupling, it 
follows that the diagonal component of a r is just 9/ir 
times the fine structure constant (which is e 2 c/x /2/i in 
SI units). 



B. Theory of orbital magnetoelectric coupling 

The purely electronic orbital magnetoelectric coupling 
OLij can be written in terms of three gauge-invariant con- 
tributions 
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(4) 



<5i,a ua is the above-mentioned (isotropic) 



where af^ _ Uyl 
CSOMP, while and a 1 ^ arc two additional contribu- 
tions. The isotropic part of the OMP tensor has contribu- 
tions from the two a terms as well as from the CSOMP 
term. The three contributions to the OMP can com- 
pactly be expressed as 
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(5) 
(6) 



a i 3 =V£jkibn I d k{d k u n \ l \D i u m \ t ){u m \<\{diH\c)\u n k), 

(7) 

where the notations are defined as follows. An implied 
sum notation applies to repeated Cartesian (ijkl) and 
band (mn) indices, corresponding to a trace over occu- 
pied bands in the latter case (written explicitly as 'tr'). A 
common prefactor 77 = — e/h(2Tr) 3 appears in each equa- 
tion, with e > being the magnitude of the electron 
charge. The Berry connection 

Annkj = (Umk\idj\Unk) (8) 

is defined in terms of the cell-periodic Bloch functions 

Kk) =e- jk - r |7/W), (9) 

which are the eigenvectors of Hy = e~ T/ He T , where 
H is the bulk periodic Hamiltonian of the crystal at zero 
electric and magnetic field, dj and Dj are the partial 
derivatives with respect to the j-th component of the 
wavevector k and the electric field S respectively. Finally, 
the tilde indicates a covariant derivative, dj = Qkdj and 
Dj = QkDj, where = 1 — |7i n k)( u nk| (sum implied 
over n). Additional screening contributions to a\p and 
<5 1( p that occur in the context of self-consistent field cal- 
culations, not given here, can be found in Ref. 11. 

As in the case of electronic polarization, one needs to 
be careful about relating the above bulk expressions to 
experimentally measurable physical quantities, since ar- 
bitrary surface modifications can contribute to the ef- 
fective measurable OMP. The relationship between the 
OMP and experimentally measurable responses are ex- 
plained in more detail in the next section. 



C. Relation between bulk and surface properties 

In order to discuss the relationship between bulk and 
surface quantities in connection with the OMP, it is in- 
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structive first to review the corresponding connections in 
the theory of electric polarization. 

1. Electric polarization and surface charge 

We first review the relationship between the bulk elec- 
tric polarization, as obtained from the crystal bandstruc- 
ture according to the Berry-phase theory, 20,21 and a mea- 
surable quantity which is the macroscopic dipolc moment 
of a finite sample cut from this crystal. Given the set of 
valence Bloch wavefunctions |-0nk} of an insulating crys- 
tal, one can readily calculate the electronic contribution 
to the polarization as the integral 

Vi = -T^E / d3k KkR, K lk ) (io) 

over the Brillouin zone (BZ). Gauge changes (|u„k) — > 
e~ l PW\u n k)) can change the value of this integral only 
by Re/f2, where R is a lattice vector and fi is the unit 
cell volume. The value of this integral is therefore only 
well-defined modulo TLe/fl. In what follows we assume 
that a definite choice of gauge has been made so that a 
definite value of V has been established. We now analyze 
how, and under what circumstances, one can relate this 
T to the (experimentally measurable) dipolc moment d 
of an arbitrarily faceted finite sample of this crystal. 

At each local region on the surface of this finite sample, 
assuming a perfect surface preparation (defect-free with 
ideal periodicity), we can relate "P to the surface charge 
density a at that same point via 21 

a = (v+ ^R) -h + A. (11) 

Here n is the surface normal unit vector, R is a lattice 
vector, and A is an additional contribution present only 
for metallic surfaces. The term involving R, which cor- 
responds to an integer number of electrons per surface 
unit cell, is required because, for a given surface n, it 
may be possible to prepare the surface in different ways 
(e.g., by adding or subtracting a layer of ions, or by fill- 
ing or emptying a surface band) such that the surface 
charge per cell changes by a quantum. Thus, R is in 
general a surface-dependent quantity in Eq. (11). If the 
surface patch under consideration is not insulating, then 
A is a term which measures the contribution of the par- 
tially occupied surface bands to the surface charge, and 
is proportional to the area fraction of occupied band in 
k space. (In the case of an insulator with non-zero first 
Chcrn number, this fraction has to be calculated with 
special care, 22 but we shall not consider this case in what 
follows.) 

Now, let us consider the special case that all surfaces 
are insulating (A = 0) and that the surface charges of all 
surface patches are consistent with a single vector value 
of R ("global consistency"). Under these circumstances, 
the macroscopic dipole moment d of the crystallite is 



given by 

d = v(7>+^R), (12) 

which can be obtained trivially by integrating Eq. (11). 
Here V is the volume of entire finite sample. As could 
be anticipated, d/V has a component depending only 
on the bulk wavefunctions and our gauge choice, and an 
additional component eR/ft reflecting the preparation of 
the surfaces. 



2. OMP and surface anomalous Hall conductivity 

We now discuss a corresponding set of relationships be- 
tween the bulk-calculated OMP and the surface anoma- 
lous Hall conductivity. 

Using Eqs. (5), (6), and (7) one can calculate the tensor 
a. from the knowledge of bulk Hamiltonian of an insulat- 
ing crystal. Analogously as in the case of polarization, 
one can again show that a gauge change 23 must either 
leave a invariant or change it by a quantum m(e 2 / 'h)l, 
where m is an integer and I is the unit matrix. More 
precisely, this gauge transformation will only affect the 
CSOMP component a cs of the OMP, since the other two 
contributions 5 LC and 5 IC are fully gauge-invariant (see 
Ref. 11 for details). 

We now imagine cutting a finite crystallite from this 
infinite crystal, and we wish to relate a. to its physically 
observable linear magnetoelectric coupling /3, defined for 
a finite sample by 

_dd L _d j i L 

l lJ ~ OB, - d£, ' 1 ' 

where di is the dipole moment of the finite sample and fij 
is its magnetic dipole moment. We want to discuss this 
relationship in a way that is analogous to that between 
the bulk T and sample dipolc moment d in Sec. II C 1. 

As follows from Eq. (1), the application of an electric 
field £j to the insulating crystal induces the magnetiza- 
tion 

M k = a jk £j, (14) 

where a is given by Eq. (4) and is only determined mod- 
ulo the quantum m(e 2 / h)l. Having a homogeneous Mk 
inside the sample and Mk = outside is equivalent to 
having a surface current Ki equal to 

Ki = €ikiM k ni, (15) 

where n\ is the surface unit normal. By eliminating Mk 
from these equations, we see that having a magnetoelec- 
tric tensor a is equivalent to having a surface anomalous 
Hall conductivity <7^ H = etkiajkrii- If the surface patch 
in question is insulating, then its anomalous Hall con- 
ductivity should just be given, modulo m(e 2 //i)I, by this 
equation. If instead the surface patch is metallic, then 



an additional surface contribution A, 3 - should be present, 
leading to the relation 



AH 



e»M I ce jk + m—Sjk ) ni + A 



(16) 



This equation is in precise analogy to Eq. (11) relating 
the polarization to the surface charge. Here Ay may in 
general contain dissipative contributions, but in the dirty 
limit it will be dominated by the intrinsic surface contri- 
bution that can be calculated as a 2D BZ integral of the 
Berry curvature of the occupied surface states. 24 The in- 
teger quantum m appearing in Eq. (16) corresponds to 
the theoretical possibility that the surface preparation 
can be changed in such a way that a surface band having 
a nonzero Chern number may become occupied. For ex- 
ample, this could be done in principle by constructing a 
2D quantum anomalous Hall layer (as described, e.g., by 
the Haldanc model 25 ), straining it to be commensurate 
with the surface, and adiabatically turning on hopping 
matrix elements to "stitch it" onto the surface. 

In the special case that all surface patches are insulat- 
ing (Ay- = 0), and all surface patches have an anomalous 
Hall conductivity given by Eq. (16) with the same value 
of m ( "global consistency" ), we can relate the experimen- 
tally measurable magnetoelectric response (3 of the finite 
crystallite to the bulk-calculated ct via 



(3 = V ( a + m—I 



(17) 



which follows by integrating Eq. (16) over all surfaces. 
This equation is in close analogy to Eq. (12) for the case 
of electric polarization. In particular, we see that f3/V 
has a component a depending only on the bulk wave- 
functions and our gauge choice, and an additional com- 
ponent that is an integer multiple of (e 2 /h)I, reflecting 
the preparation of the surfaces. 

As will be discussed in the next section, time-reversal 
symmetry imposes additional constraints on a, and some 
care is needed in the interpretation of Eq. (17) for the 
case of Z2 topological insulators. 



D. Motivation and relationship to strong Z2 
topological insulators 

In this Section, we give arguments to motivate our 
hope that in certain materials the CSOMP might be 
on the order of, or even much larger than, the total 
magnetoelectric coupling in typical known magnetoelec- 
tric materials. For simplicity, we focus henceforth only 
on the CSOMP part of the total OMP response, even 
though there are additional contributions coming from 
5 LC and 5 IC . Thus, from now on, the quantity 9 mea- 
sures the strength of the CSOMP through the relation 
a cs = 9e 2 /2Tih. 
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FIG. 1. Identical samples cut from a strong Z2 topological 
insulator, but with two different surface preparations, (a) 
Time-reversal symmetry is preserved at vacuum-terminated 
surfaces; the net magnetoelectric coupling of this sample is 
zero, (b) Time-reversal symmetry is broken at the surface as 
a result of exchange coupling to an insulating ferromagnetic 
adlayer; if this opens a gap in the surface-state spectrum, 
the entire sample will behave as if it has a magnetoelectric 
coupling of exactly 8 = tv. 



1. Time-reversal symmetry constraints on 9 

Let us analyze the allowed values of 9 for an infi- 
nite bulk insulating system that respects time-reversal 
(T) symmetry. Since T flips the sign of the magnetic 
field, it will also reverse the sign of 9. As mentioned 
earlier in Sec. II C 2, however, the value of 9 can be 
changed by 27r under a gauge transformation. Therefore 
one concludes 12 ^ 13 that the allowed values of 9 consis- 
tent with T symmetry are mod 2-7T and it mod 27T, and 
that these two cases provide a topological classification of 
all T-invariant insulators. Indeed, this classification has 
been shown 12,13 to be identical to the one based on the 
Z2 index, with Z2-odd or "strong topological" insulators 
having 9 — n, while Z2-even or "normal" insulators have 
9 = 0, even though the Z2 index is most often introduced 



in a different context. 6 (Incidentally, a 



in 



both cases since these terms are fully gauge-independent, 
unlike the CSOMP term which can be changed by 27r.) 

Consider now a finite sample of a normal (Z2-even) T- 
symmetric insulator (9 = in the bulk) with insulating 
surfaces (Ay = 0) prepared in a way that the integer m 
is nonzero, and the same on every surface. From Eq. (17) 
we conclude that this sample will have a non-zero mag- 
netoelectric response, f3, proportional to m. Obviously a 
sample that has T symmetry both in the bulk and on the 
surface must have (3 = 0, and therefore we conclude that 
this system needs to have broken T reversal symmetry 
at the surface. As mentioned earlier, one could, at least 
formally, prepare such a surface by starting from the one 
that has m = and then absorbing to each surface a 
layer of anomalous Hall insulator 25 with Chern index m. 
Such a procedure will keep the surfaces insulating but it 
will necessarily break the T-reversal symmetry. 

Next we analyze the case of a strong Z2 topological 
insulator having 9 — 7r, or cquivalently, a = a cs = 
(e 2 /2h)I. We first consider a sample of such a system 
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that has T symmetry conserved at its surfaces, as in 
Fig. 1(a). Again, since the entire sample is T-symmetric, 
its experimentally measurable magnctoelectric coupling 
tensor (3 clearly has to vanish. Using Eq. (16) and the 
fact that m can take on only integer, and not half-integer, 
values, we conclude that the only way to make the re- 
sponse of the entire sample vanish is to have be non- 
zero. This requires that the surfaces of such a system 
must be metallic. Moreover, since the contribution A^- 
of the metallic surface band to the surface anomalous Hall 
conductivity is just given by the Berry phase around the 
Fermi loop, 24 the needed cancellation requires this Berry 
phase to be exactly ±7r. All this is in precise accord 
with the known properties of Z2-odd insulators and their 
topologically protected surface states. 26 

The Kramers degeneracy at the Dirac cone in the sur- 
face bandstructure can be removed by the application of 
a T-breaking perturbation to the surface. In principle, 
this could be accomplished, for example, by applying a 
local magnetic field to the surface or by interfacing the 
surface to an insulating magnetic overlayer. In the latter 
case, the interatomic exchange couplings provide a kind 
of effective magnetic field acting on the surface layer of 
the topological insulator. If the local Fermi level resides 
in the gap opened by field, then the surface becomes in- 
sulating. If the field can be consistently oriented (see 
Ref. 12) on each patch of the surface, either along or op- 
posite the direction of surface normal vector n (as shown 
in Fig. 1(b)), then the entire surface becomes insulating. 
It is important that the field is applied consistently in 
the same direction with respect to n, since conducting 
channels will otherwise appear at domain boundaries. 26 

If all of these requirements are met, the surface contri- 
bution Aij to (3 vanishes, so that (3 = Vol with a given 
only by bulk value of 9 = tt (assuming m = for simplic- 
ity). Therefore such a sample of a strong Z2 topological 
insulator would behave as if the entire sample has exactly 
half a quantum of magnetoelectric coupling (9 = n), even 
though its bulk is time-reversal symmetric! 



2. Prospects for large-8 materials 

Recently surface-sensitive ARPES measurements have 
experimentally confirmed that several compounds, 15-17 
including Bii-^Sba;, Bi2Se3, Bi2Te3 and Sb2Te3, do in- 
deed behave as strong Z2 topological insulators. There- 
fore their bulk wavefunctions must be characterized by 
8 = tt. Up to now, the corresponding magnetoelectric 
response has not been measured experimentally, in part 
because of the difficulties in obtaining truly insulating be- 
havior in the bulk, as well as the need to gap the surfaces 
by putting them in contact with magnetic overlayers as 
described earlier. 

We believe that a more promising approach to observ- 
ing a large CSOMP (i.e., 9 comparable to tt) is to consider 
an insulator that has neither T nor spatial inversion sym- 
metry. In this case the Z2 classification does not apply, 
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FIG. 2. (Color online) Schematic view of the allowable values 
of 8 in different parts of the two-parameter space of some un- 
specified model Hamiltonian. Horizontal axis corresponds to 
the perturbation that preserves at least one of the symmetries 
that render 9 to be or tt (see Sec. HE). Vertical axis pa- 
rameterizes a perturbation that breaks those symmetries and 
allows 9 to be arbitrary. See text for the details. 



and the surface can be gapped without any need to apply 
a T-brcaking perturbation. (A more precise statement of 
the symmetry considerations will be given in Sec. HE.) 
The sample can then display a bulk magnetoelectric cou- 
pling of the simple form (3 = Vol. We note that an or- 
bital magnetoelectric coupling of 9 ~ tt (i.e., a r ~ 1/137) 
would correspond to a EH ~ 24.3 ps/m, a value that is 
significantly larger than the observed coupling in Cr 2 3 , 
one of the best-studied magnetoelectric materials. For 
comparison, the reported experimental values for a'p in 
O2O3, which are presumably dominated by spin-lattice 
coupling, range between 0.7 and 1.6ps/m at 4.2 K. 27,28 

Of course, in order to have a good chance of finding 
a material with a large 9, it may be advisable to look 
for materials with some of the same characteristics as 
the known Z2-odd insulators, of which the most impor- 
tant is probably the presence of heavy atoms with strong 
spin-orbit coupling. We see no strong reason why such a 
search might not reveal a material having a large OMP 
in the above sense. 

To illustrate the kind of a search we have in mind, con- 
sider some model Hamiltonian that depends on two pa- 
rameters, one that preserves either the T or spatial inver- 
sion symmetry (or both), and another that that breaks 
symmetry such that 9 takes a generic value. The possi- 
ble behavior of such a model is sketched in Fig. 2, where 
these two parameters are plotted along the horizontal 
and vertical axes respectively. The figure also indicates 
the generic value of 9 in each region of parameter space. 
Along the horizontal axis, where the extra symmetry is 
present, three regions are indicated. The black dot in- 
dicates a point of gap closure forming the boundary be- 
tween a normal T-symmctric insulator regime on the left 
(9 = 0) and a strong Z2 topological insulator regime on 
the right (9 = tt). If the system is carried along the hori- 
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TABLE I. Magnetic point groups for which a non-zero 
CSOMP is allowed by symmetry. Notation follows Ref. 29. 
Point groups in bold allow only for a purely isotropic magne- 
toelectric tensor. 
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III. METHODS 

In this section we present our methods for calculat- 
ing the CSOMP in the framework of density-functional 
theory, and analyze in more detail its mathematical prop- 
erties and the formal similarities to the formulas used to 
calculate electric polarization and anomalous Hall con- 
ductivity. 



A. Review of Berry formalism 



zontal axis, 9 must be either or it except at the critical 
point, and it must therefore jump discontinuously when 
passing through this point of metallic behavior. On the 
other hand, if we now imagine passing from the Z2-odd 
to the Z2-even phase along the dashed curve in Fig. 2, 
9 can vary smoothly and continuously from ir to with- 
out any gap closure anywhere along the path. If we can 
identify a material lying near, but not at, the right end of 
this dashed path, it could be the kind of large-6* material 
we seek. 

Thus, our ultimate goal is to use first-principles cal- 
culations to search for a large 9, not in a topological in- 
sulator, but in an "ordinary" (but presumably strongly 
spin-orbit coupled) insulating magnetic material. While 
our work has yet to result in the identification of a large- 
9 material of this kind, it represents a first step in the 
desired direction. 



E. General symmetry considerations 

Recall that 9 is a pseudoscalar that changes sign un- 
der time-reversal and spatial-inversion symmetries (since 
B changes sign under T while £ changes sign under in- 
version). On the other hand, 9 is invariant under any 
translation or proper rotation of a crystal. Therefore if 
the magnetic point group of a crystal contains an element 
that involves T, possibly combined with a proper rota- 
tion, the value of 9 is constrained to be or n (modulo 
2tt) as discussed earlier. The same happens if the mag- 
netic point group contains inversion symmetry or any 
other improper rotation. 

All 32 of the 122 magnetic point groups that do not 
contain such symmetry elements, and which therefore 
allow for an arbitrary value of 9, are listed in Table I. 
(The bold entries in the table are those magnetic groups 
for which the tensor a must be isotropic, i.e., a constant 
times the identity matrix; the same magnetic groups were 
also analyzed in Ref. 18). Clearly we can constrain our 
search for interesting materials to the cases listed in the 
Table. 



Assume we are given the Bloch wavefunctions \ipnk) = 
e lkr |M„k) as a function of wavevector k in the d- 
dimensional BZ (d = 1, 2, or 3) for an insulator having 
valence bands indexed by n £ {1, . . . , iV}. We work with 
the cell-periodic Bloch functions u n k(r) = e~ lk ' r V>nk(r) 
and allow them to be mixed at each k point by an arbi- 
trary fc-dependent unitary matrix 

Kik) ->■ \u m k)U mn k (18) 

(sum on m implied). After this gauge transformation the 
wavefunctions are no longer eigenfunctions of the Hamil- 
tonian, but they span the same A^-dimensional subset of 
the Hilbert space as the true eigenfunctions. For any 
given choice of gauge, we define the Berry connection 

d 

which is a fc-dependent N x N x d matrix that measures, 
at each k point, the infinitesimal phase difference be- 
tween the m-th and n-th wavefunctions associated with 
neighboring points along Cartesian direction j in k space. 
This object was already briefly introduced in Eq. (8). 

In the context of electronic structure calculations, we 
can now list three material properties that can be eval- 
uated knowing only the Berry connection: the electric 
polarization, the intrinsic anomalous Hall conductivity, 
and the CSOMP. 

The electric polarization V already appears in dimen- 
sion d = 1 and it can be evaluated as an integral of the 
Berry connection over the one-dimensional BZ as 20 

V = ~ [ dktvAk, (20) 

JBZ 

where the trace is performed over the band indices of the 
Berry connection, as in Eq. (10). The integrand is also 
referred to as the Chcrn-Simons 1-form, and its integral 
over the BZ is well known to be defined only modulo 
27r. Any periodic adiabatic evolution of the Hamiltonian 
%(A) whose first Chern number in (k, A) space is non-zero 
will change the integral above by a multiple of 27r. 20 

Unlike onc-dimcnsional systems, crystals in d = 2 can 
have an anomalous Hall conductivity. For a metal, the in- 
trinsic contribution from a band crossing the Fermi level 
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FIG. 3. Graphical interpretation of Eqs. (20) (a), (21) (b) 
and (22) (c) in the case of one occupied electron band and 
for cubic crystal symmetry, for simplicity. See text for more 
detail. 



can be evaluated as a line integral 



24,30 



.AH 



e 2 1 

= — — cf> dk ■ Ak 



h 2tt 



(21) 



FL 



over the Fermi loop. Fully-filled deeper bands can also 
make a quantized contribution given by a similar inte- 
gral, but around the entire BZ; this is the only contribu- 
tion in the case of a quantum anomalous Hall insulator. 25 
(In both cases, the gauge choice on the boundary of the 
region should be consistent with a continuous, but not 
necessarily fc-periodic, gauge in its interior; alternatively, 
each expression can be converted to an area integral of a 
Berry curvature to resolve any uncertainty about branch 
choice. See Rcf. 31 for more details.) 

Finally, unlike one- or two-dimensional systems, three- 
dimensional systems can have an isotropic magnetoelec- 
tric coupling. The CSOMP can be evaluated in d = 3 as 
a BZ integration of a quantity involving the Berry con- 
nection: 



47T 



d keijktT 



BZ 



AiOjAk g AiAjAk 



(22) 



The integrand in this expression is known as the Chcrn- 
Simons 3-form, and its integral over the entire BZ is again 
ill-defined modulo 2-7T, since any periodic adiabatic evolu- 
tion of the Hamiltonian H(X) whose second Chern num- 
ber in (k, A) space is non-zero will change 9 by an integer 
multiple of 2ir. 12,13 

The sketches in Fig. 3 compare the geometrical char- 
acters of the operations needed to evaluate Eqs. (20-22) 
in practice. We consider the case of one occupied elec- 
tron band for simplicity. The polarization of Eq. (20) is 
calculated by a line integral; on a discrete fc-mesh, the 
integral of the Berry connection A over each of line seg- 
ment, as in Fig. 3(a), is converted to a discretized form 
(see Eq. (23)). Similarly, in two dimensions the anoma- 
lous Hall conductivity of Eq. (21) can be calculated as 
suggested in Fig. 3(b) by dividing the occupied part of 
the BZ into small square segments and then integrating 
A around each square. (Equivalently, one can integrate 
A along the Fermi loop. 31 ) In three dimensions, Fig. 3(c), 
Eq. (22) can be evaluated by dividing the BZ into small 
cubes. In each, one needs to multiply the integral of 
A along one of the Cartesian directions (as in Eq. (20)) 



with the integral of Berry connection in the square or- 
thogonal to that direction (as in Eq. (21)), followed by a 
symmctrization over the three Cartesian directions. 



B. Numerical evaluation of 6 

In electronic-structure calculations, the cell-periodic 
wavefunctions |u„k) are typically calculated on a uniform 
fc-space grid with no special gauge choice; in general, one 
should assume that the phases have been randomly as- 
signed. Nevertheless, it is straightforward to construct a 
gauge-invariant polarization formula that is immune to 
this kind of scrambling of the gauge. 32 In one dimension 
with kj for j G {1, . . . , M} (where km is the periodic im- 
age of point fci), the electronic polarization is calculated 
as 



P = — Imlndet [M klk2 M k2k3 ...M kM _ lkM ] 
where the overlap matrix M kk > is defined as 



\M, 



kk> 



{V"mk \ ^nk' ) • 



(23) 



(24) 



The reason for using Eq. (23) is that the determinant 
of the matrix M klk2 M k2ka ...M kM ^ lkM is gauge-invariant 
under any transformation in the form of Eq. (18). Ad- 
ditionally, the implementation of Eq. (23) is numerically 
stable even when there are band crossings. A similar 
gauge-invariant discretization can also be used to calcu- 
late the anomalous Hall conductivity er AH . 31 

Unfortunately, except in the single-band ("Abelian") 
case, we are unaware of any corresponding gauge- 
invariant discretized formula for the integral of the 
Chern-Simons 3-form. As a result, we have no prescrip- 
tion for computing the CSOMP that is exactly gauge- 
invariant for a given choice of k mesh. This is a se- 
rious problem. Unlike the calculation of the polariza- 
tion, which is straightforward even if the gauge is ran- 
domly scrambled at each mesh point, the calculation of 
the CSOMP requires that we first identify a reasonably 
smooth gauge on the discrete mesh. 

The problem of finding a smooth gauge in k is essen- 
tially the same as that of finding well-localized Wannicr 
functions. For this reason, we have adopted here the 
approach of first constructing a Wannicr representation 
for the valence bands, and then using it to compute the 
CSOMP. In fact, starting from Eq. (22), we derive an ex- 
pression that allows us to compute 9 directly in the Wan- 
nier representation. Once we have well-localized Wan- 
nier functions, this guarantees smoothness of the gauge 
and avoids problems with band crossings. Admittedly, 
such a formula still depends on the gauge choice, mean- 
ing that different choices of Wannicr functions will lead 
to slightly different results. However, this difference will 
vanish as one increases the density of the fc-point mesh, 
since in the continuum limit the /c-space expression for 
9 is gauge-invariant (modulo 2tt). More precisely, we ex- 
pect the calculation of 9 to converge once the inverse of 
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the /c-point mesh spacing becomes much larger than the 
spread of the Wannier functions. 

Therefore, we adopt the strategy of calculating 9 on 
k meshes of different density, and extrapolating 9 to the 
limit of an infinitely dense mesh. Furthermore, we con- 
struct maximally-localized Wannier functions (MLWF) 
following Ref. 32, expecting this to give relatively rapid 
convergence as a function of the k mesh density. 

Recall that the Wannier function associated with (gen- 
eralized) band index n in unit cell R is defined in terms 
of the rotated Bloch states (18) as 



R/ 



CI 



(2tt) 3 



<Fk 



ik-(r-R) 



(25) 



In the case of MLWFs, the C/ m „k are chosen in such a way 
that the total quadratic spread of the Wannier function 
is minimized. 32 (In practice the BZ integral is replaced 
by a summation over a uniform grid of k points.) 

Using Eq. (25), one can relate the Berry-connection 
matrix A mn kj in the smooth gauge to the Wannier ma- 
trix elements of the position operator through 32 



A m nkj =^e ik ' R (Om|rj|Rn). 



(26) 



R 



Replacing each occurrence of Aj in Eq. (22) with the 
above gives, after some algebra, 



, 1 (2tt) ; 
"4tt fl 



e ilk lm 



- ^(Om\r i \'Rn){Rn\r j \Orri)R k 



R 



-^(Oi|r<|Rro)<Rm|r 3 -|Pn)(Pn|r fc |(H) 



RP 



(27) 



where the sum is implied over band (Iran) and Cartesian 
(ijk) indices. 

To obtain a more symmetric form, we introduce a mod- 
ified position-operator matrix element between WFs de- 
fined as 



(Rm\n\Pn) = (Rm|r.;|Pn) (1 - 5 mn S RF ) 
and a notation for the Wannier center 

r ni = (On | Ti| On). 
Then Eq. (27) becomes 

„ 1 W 



4tt n 



"Eii'fc X 



(28) 
(29) 

(30) 



Illl 



(Om | h | Rn) (Rn | fj \ Om) (Rk + r nk 



Tmkj 



R 



^2 -{Ol\fi\Km)(Rm\f j \Pn)(Pn\r k \Ol) 



RP 



(31) 



We find this form more convenient because it separates 
the contributions of diagonal and off-diagonal elements 
of position operators. 33 (It is also manifestly invariant to 



the reassignment of a Wannier function to a neighbor- 
ing cell.) The validity of Eqs. (27) and (31) has been 
tested numerically by comparing with the evaluation of 
Eq. (22) for the case of a tight-binding model introduced 
in Ref. 11. The evaluated expressions agreed to numer- 
ical accuracy after extrapolation to the infinitely dense 
mesh. These expressions can also be shown to be gauge- 
invariant by working directly within the Wannier repre- 
sentation. 



C. Computational details 

Calculations of the electronic ground state and 
of structural relaxations were performed using the 
Quantum-ESPRESSO package, 34 and the Wannier90 
code 35 was used for constructing maximally localized 
Wannier functions. We used radial-grid discretizcd 
HGH 36 norm-conserving pseudopotcntials. Calculations 
were performed in the noncollincar spin framework, in- 
cluding spin-orbit effects as incorporated in the pseu- 
dopotentials. In all calculations we used the Perdew- 
Wang 37 LDA energy functional. The pseudopotentials 
used for Cr, Fe and Gd contain semi-core states, while 
the ones for Al, Bi, Se and O do not. 

The self-consistent calculations on Cr203 were per- 
formed on a 4 x 4 x 4 Monkhorst-Pack 38 grid in k space. 
Non-self-consistent calculations for the Wannier-function 
construction were performed on /c-space grids containing 
the origin and ranging in size from 6x6x6 to 12x12x12. 
The plane-wave energy cutoff was chosen to be 150 Ry. 

In the case of Bi2Se3, the self-consistent calculations 
were performed on a 6 x 6 x 6 grid with energy cutoff of 
60 Ry, while the non-selfconsistent calculation was done 
on grids between 6x6x6 and 11 x 11 x 11. 

The position-operator matrix elements (Om|rj|Rn) 
needed to evaluate Eq. (31) were calculated in k space 
by inverting the Fourier sum in Eq. (26) over the non- 
self-consistent /c-point mesh, and then approximating the 
k derivative in Eq. (19) by finite differences on that mesh, 
as detailed in Ref. 30. 



IV. RESULTS AND DISCUSSION 

A. Conventional magnetoelectrics 

In this section we present the results of our first- 
principles electronic-structure calculations of 9. We be- 
gin with conventional magnetoelectrics, i.e., materials 
that are already experimentally known to have a non- 
zero magnetoelectric tensor. Some of these materials do 
not allow all diagonal components of the magnetoelec- 
tric tensor to be non-zero. We omit those materials from 
our analysis here, since we are interested in calculating 
the CSOMP part of the magnetoelectric coupling, which 
would vanish in such cases. We first present our results 
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FIG. 4. (Color online) (a) Rhombohedral unit cell of Cr2 03. 
Magnetic moments on Cr atoms are indicated by red arrows 
and oxygen octahedra are drawn around each Cr atom, (b) 
Schematic of hexagonal unit cell of Bi2Se3 with imposed local 
Zeeman field on Bi atoms. Induced magnetic moments are 
shown by red arrows. Thick blue lines indicate Se layers; 
letters (ABC) indicate stacking sequence. In both panels, the 
vertical line indicates the 3-fold rhombohedral axis, and the 
cross designates a 2-fold rotation axis orthogonal to the plane 
of the figure (also a center of inversion coupled with time 
reversal) . 



on Ci'203 in some detail, and then briefly discuss our 
results for BiFe0 3 and GdA10 3 . 




0.0 



0.00 



0.05 



0.10 0.15 

Ak (A" 1 ) 



0.20 



0.25 



FIG. 5. (Color online) Calculated value of 9 in Cr2 03 for 
varying densities of fc-space grids, where Afc is the nearest- 
neighbor distance on the grid. Top axis specifies the size of the 
corresponding uniform Monkhorst-Pack grid. Line indicates 
a quadratic extrapolation of 9 to the infinitely dense k mesh. 



line indicates the second-order polynomial extrapolation 
to an infinitely dense mesh. The extrapolated value of 9 
is 1.3 x 10~ 3 , which is a small fraction of the quantum of 



OMP 9 = 2tt and corresponds to a*j£ = a™ = a™ = 
0.01 ps/m. The positive sign of 9 pertains to the pattern 
of Cr magnetic moments shown in Fig. 4(a); reversal of 
all magnetic moments would flip the sign of 9. 

In order to compare this value of the magnetoelcctric 
coupling with experimental values and other theoretical 
calculations, we somewhat arbitrarily define 



Calculation of 9 in O2 O3 



Q 



off 



n. 



(32) 



We first fully relax the structure in the R3c space group 
and obtain the Wyckoff position to be x = 0.1575 for Cr 
atoms (4c orbit) and x = —0.0690 for O (6c orbit). The 
length of the rhombohedral lattice vector is a = 5.3221A 
while the rhombohedral angle is 53.01°. The Cr atoms 
have magnetic moments pointing along the rhombohedral 
axis as illustrated in Fig. 4(a) in an antiferromagnetic ar- 
rangement. The value of the magnetic moment is 2.0 /j,b 
per Cr atom and the electronic gap is 1.3 eV, which agrees 
well with previous LDA+U calculations 39,40 in the limit 
where the on-site Coulomb parameter U is set to zero. 

Neglecting for a moment the magnetic spins on the Cr 
sites, the space-group generators are a three-fold rota- 
tion, a two-fold rotation, and an inversion symmetry as 
indicated in Fig. 4(a). Its point group is therefore 3m. 
If we now include the spins on the Cr atoms in the anal- 
ysis, we find that the three-fold and two-fold rotations 
remain, while the inversion becomes a symmetry only 
when combined with time-reversal. Therefore the mag- 
netic point group of Cr203 is 3'm'. 41 This magnetic point 
group allows 9 to be different from or 7r, as discussed 
in Sec. HE. 

Figure 5 shows the calculated values of 9 using Eq. (31) 
for &2O3 with A:-space meshes of various densities. The 



The value of a cS obtained from the results of Ref. 8 
is 0.23 ps/m for the purely electronic part of the spin- 
mediated component. Therefore, our calculated CSOMP 
contribution in C^Os amounts to only 4% of this elec- 
tronic spin component. The ionic component of the 
spin response calculated by the same authors results in 
a cS = 0.74 ps/m, while the one calculated in Ref. 7 is 
about 2.6 times smaller, 0.29 ps/m. (In both of these 
calculations, a zz is zero.) Finally, experimental mea- 
surements of the magnetoelectric tensor in Cr203 at 
4.2 K vary between a = 0.55 ps/m and 1.17ps/m (see 
Refs. 27 and 28 respectively). 

Clearly, our computed CSOMP contribution for Cr 2 03 
is negligible, being two orders of magnitude smaller than 
the dominant lattice-mediated spin contribution. This is 
probably not surprising, since the spin-orbit coupling is 
relatively weak in this material. Given that it is weak, 
we can guess that that magnitude of the CSOMP should 
be linear in the strength of the spin-orbit interaction in 
Cr203. Our calculations allow us to check this by vary- 
ing the spin-orbit interaction strength Aso between (no 
spin orbit) and 1 (full spin-orbit interaction). As shown 
in Fig. 6, if we calculate 9 for various intermediate val- 
ues of Aso, we see that the CSOMP docs indeed depend 
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FIG. 6. (Color online) Calculated 8 in Cr203 as a function 
of spin-orbit coupling strength, scaled such that Aso = 1 cor- 
responds to the full spin-orbit coupling strength and 9o = 
6>(Aso = !)■ 
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FIG. 7. (Color online) Calculated value of 6 in Bi2Se3 for 
varying densities of fc-space grids, where Afc is the nearest- 
neighbor distance on the grid. Top axis specifies the size of the 
corresponding uniform Monkhorst-Pack grid. Line indicates 
a quadratic extrapolation of 9 to the infinitely dense k mesh. 



roughly linearly on Ago- 



2. Other conventional magnetoelectrics 

We have also carried out calculations of 9 in BiFcC>3 
and GdAlC>3 , but with a smaller number of fc-point grids 
than in the case of Cr2C>3. Therefore, our results are 
less accurate, but should still give a correct ordcr-of- 
magnitude estimate of 6. 

For BiFeOa we perform the calculation in the 10-atom 
antifcrromagnetic unit cell (the long- wavelength spin spi- 
ral was suppressed). We obtain an electronic band gap 
of 0.95 eV with magnetic moments of 3.5 /^b on each Fe 
atom, and with a net magnetization of 0.1 [1b per 10-atom 
primitive unit cell due to the canting of the Fe magnetic 
moments. Extrapolating 9 to an infinitely dense mesh 
using just 6x6x6 and 8x8x8 fc-point meshes, we ob- 
tain = 0.9 x 10~ 4 . In the case of GdAlOs we calculate 
the electronic band gap to be 5.0 eV and the Gd magnetic 
moment to be 6.7 /xb- We obtain a value of 9 = 1.1 x 10~ 4 
after extrapolating calculations using 4x4x4 and 6x6x6 
fc-space meshes. Thus, it is clear that the CSOMP is very 
small in both materials. 



B. Strong Z2 topological insulators 

We now investigate the CSOMP in the case of Bi2Se3, 
which is known experimentally 17 and theoretically 42 to 
belong to the class of strong Z2 topological insulators. 
In the absence of broken T symmetry, such a material 
should have a 9 of exactly it (modulo 2ir). We first con- 
firm this numerically. Then, in Sec. IV C, we also study 
what happens when T is broken artificially by inducing 
antiferromagnetic order on the Bi atoms and tracking the 
resulting variation of 9. 

Bi2Sc3 is known to belong to space group R3m, with 



Bi at a 2c site and Se at the high-symmetry la site as 
well as at a 2c site. In our calculations we find that 
the Wyckoff parameters for Bi and Se are x = 0.4013 
and 0.2085 respectively. We also find the length of the 
rhombohedral lattice vector to be a = 9.5677 A and the 
rhombohedral angle to be only 24.77°. The electronic 
gap is calculated to be 0.4 eV. 

The generators of the R3m space group are again three- 
fold and two-fold rotations and inversion (point group 
3m). Since the system is nonmagnetic, the magnetic 
space group also contains the T symmetry operator, and 
its magnetic point group is 3ml'. According to the anal- 
ysis given in Sec. HE, it is clear that 9 must therefore be 
zero or 7r (modulo 2-7r). 

Since we know that Bi2Sc3 is a strong Z2 topologi- 
cal insulator, we expect that 9 should be equal to ti 
(modulo 2ir). However, special care needs to be taken 
in order to evaluate 9 in such a case, because the choice 
of a smooth gauge becomes problematic. Specifically, 
it is known that the Z2 topology presents an obstruc- 
tion to the construction of a Wannier representation (or 
equivalently, a smooth gauge in k space) that respects T 
symmetry. 43,44 Therefore, during the maximal localiza- 
tion procedure, one needs to choose trial Wannier func- 
tions that do not take the form of Kramers pairs, thereby 
explicitly breaking the T symmetry. 45 (It is important to 
note that this choice of Wannier functions does not bias 
our calculation towards having 9 = tt, since the same 
starting choice of T-symmetry-broken Wannier functions 
for a normal T-symmetric insulator would result in 9 = 
up to the numerical accuracy of the calculation.) 

Our results for 9 in Bi2Se3 are given in Fig. 7 for var- 
ious densities of k meshes, ranging from 6 x 6 x 6 to 
llxllxll. A quadratic polynomial extrapolation to 
the infinitely dense mesh limit gives 9 — 1.0 77r. This is in 
reasonable agreement with the expected value of 9 = n, 
given the uncertainties in the extrapolation. (Of course, 
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FIG. 8. Calculated value of (vertical axis) and induced 
magnetic moment on the Bi atom (horizontal axis) for Bi2Se3 
with artificially applied staggered Zeeman field on Bi atoms, 
as described in the text. 9° is the value of CSOMP when 
magnetic field is not present. 



if we make a time-reversed choice of starting Wannier 
functions, we obtain 8 = — 1.077T, which is consistent, 
within the errors, with 8 = —n and modulo 2ir to 8 = n.) 
Clearly the convergence with respect to mesh density is 
somewhat slow, making a precise extrapolation difficult. 
The reasons for this, and some possible paths to improve- 
ment, will be discussed in Sec. V. 



C. Z2-derived nontopological insulators with 
broken symmetries 

Even though 8 = tt in Bi2Sc3, a finite sample with T 
symmetry preserved everywhere, including at the sur- 
faces, will not exhibit any magnetoelectric coupling. 
From the point of view of the discussion in Sec. II D, this 
happens because of an exact cancellation between 8 = ±tt 
contributions coming from the bulk (a) and metallic sur- 
face (A) terms in Eq. (16). However, if one breaks the 
T symmetry in the bulk (and possibly some other bulk 
symmetries, as detailed in Sec. HE), the CSOMP term 
can become allowed. 

The magnetic space group of Bi2Sc3 contains both T 
and spatial inversion symmetries. The presence of either 
by itself is enough to insure that 8 = or tt (modulo 
2tt). Now let us consider turning on, "by hand," a local 
Zeeman field on each Bi atom in the staggered arrange- 
ment shown in Fig. 4(b), i.e., with fields oriented parallel 
to the rhombohcdral axis and alternating in sign. The 
induced magnetic moments along the three-fold axis pre- 
serve both three-fold and two-fold rotation symmetries; 
both inversion and T symmetries are broken, but T taken 
together with inversion is still a symmetry. The resulting 
magnetic point group of the system is again 3'm', as it 
was for Cr203, and it does allow for a CSOMP (the same 
magnetic arrangement has also been discussed in Ref. 46 
in a different context). 



In the density functional calculation one can easily ap- 
ply a local Zeeman field to individual atoms in an ar- 
bitrary direction. 47 Using this method, we have calcu- 
lated the CSOMP in Bi2Sc3 with the pattern of local 
fields described previously and illustrated in Fig. 4(b). 
Fig. 8 presents the calculated values of 8 as a function 
of induced magnetic moment on Bi, where a positive /iBi 
corresponds to the pattern of magnetic moments indi- 
cated in Fig. 4(b). (Actually this was done by applying 
the full extrapolation procedure of Fig. 7 for one case, 
//Bi = 0.16 /Xbj and using this to scale the results cal- 
culated on the 10 x 10 x 10 grid at other /zsi-) The 
dependence of the change in CSOMP on the magnetic 
moment is linear over a wide range. One can see that 
for a relatively moderate magnetic moment of ±0.27 /^b, 
the value of 8 is changed from n to 7r±0.55. (For much 
higher local magnetic fields, Bi2Se3 becomes metallic and 
the CSOMP becomes ill-defined.) 

These results indicate that it is possible, at least in 
principle, for a magnetic material to have a large but 
unquantized value of 8, thereby providing an incentive for 
future searches for materials in which such a state arises 
spontaneously, without the need to apply perturbations 
by hand as done here. 



V. SUMMARY AND OUTLOOK 

In this manuscript, we have presented a first-principles 
method for calculating the Chern-Simons orbital magne- 
toelectric coupling in the framework of density-functional 
theory. We have also carried out calculations of this 
coupling for a few well-known magnetoelectric materi- 
als, namely Cr203, BiFeOs and GdA103. Unfortunately, 
in these materials the CSOMP contribution to the to- 
tal magnetoelectric coupling is quite small. This is not 
surprising, since in most magnetoelectric materials the 
coupling is expected to be dominated by the lattice- 
mediated response, whereas the CSOMP is a purely elec- 
tronic (frozen-ion) contribution. Moreover, the CSOMP 
is part of the orbital frozen-ion response, which is again 
expected to be smaller then the spin response, except per- 
haps in systems with very strong spin-orbit coupling, as 
discussed in Sec. I. For example, in Cr 2 03 the CSOMP 
is about 4% of the frozen-ion spin contribution to the 
magnetoelectric coupling. 

On the other hand, we have reasons to believe that 
in special cases the CSOMP contribution to the magne- 
toelectric coupling could be large compared to the total 
magnetoelectric coupling in known magnetoelectrics such 
as Cr203. After all, as already pointed out in Sec. II D 2, 
Z2 topological insulators are predicted to display a large 
magnetoelectric effect of purely orbital origin when their 
surfaces are gapped in an appropriate way. If this is so, 
why shouldn't a similar effect occur in certain T-broken 
systems? 

As a proof of concept for the existence of those special 
cases, we have considered Bi2Se3 with inversion and time- 
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reversal symmetries explicitly broken "by hand." Here 
we find that with a relatively modest induced magnetic 
moment on the Bi atoms, one can still achieve quite a 
large change in the CSOMP. 

On the computational side, there still remain several 
challenges. For example, the convergence of our calcu- 
lations of the CSOMP with respect to the fc-point mesh 
density is disappointingly slow. A direct calculation of 
9 in Bi2Se3 using a very dense mesh of 11 x 11 x 11 k 
points only manages to recover about 30% of the con- 
verged value of 9 = ir, and an extrapolation procedure 
is needed to brings us within 10% of that value. This 
clearly points to the need for methodological improve- 
ments, and we now comment briefly on some possible 
paths for future work. 

The slow convergence that we observe is related in part 
to the way in which we evaluate the position-operator 
matrix elements (Om\rj\TLn). As discussed in Ref. 30, 
the /c-space procedure we adopted (see Sec. IIIC) entails 
an error of 0(Ak 2 ). Preliminary tests on a tight-binding 
model suggest that an exponentially fast convergence of 
9 can be achieved by an alternative procedure, in which 
the WFs are first constructed on a real-space grid over 
a supcrccll (whose size scales with the fc-mesh density), 
and the position matrix elements are then evaluated di- 
rectly on that grid, as in Ref. 48. It may also be possible 
to improve the /c-space calculation by using higher-order 
finite-difference formulas that have a more rapid conver- 
gence with respect to mesh density. 

An alternative approach would be to develop a formula 
for the CSOMP that is exactly gauge invariant in the case 
of a discretized /c-space grid. Such an expression already 
exists for the case of electronic polarization, Eq. (23), but 
we are aware of no counterpart for the CSOMP. Even 
though such an approach would not necessarily provide 
much faster convergence with respect to the fc-space sam- 
pling, it would still be a significant improvement. For ex- 
ample, one would not need to construct a smooth gauge 
in k space, which is a particular problem in the case of 
Z2 insulators (or for a symmetry-broken insulator in the 
vicinity of a Z2 phase). Another use of such a formula 
would be to calculate with relative ease the Z2 index of 
any insulator, even in the cases when other methods 49-51 
cannot be applied (for example, when inversion symme- 
try is not present). 



Furthermore, a full calculation of the electronic con- 
tribution to the orbital magnetoelectric response should 
also include the remaining two contributions given in 
Eqs. (6) and (7). This calculation would also require a 
knowledge of the first derivatives of the electronic wave- 
functions with respect to electric field. While these 
derivatives are available as part of the linear-response 
capabilities of the Quantum-ESPRESSO package, 34 
some care is needed to arrive at a robust implementa- 
tion of Eqs. (6) and (7), as will be reported in a future 
communication. 

Finally, recall that our calculations have all been car- 
ried out in the context of ordinary density-functional the- 
ory. In cases where orbital currents play a role, it is 
possible that current-density functionals 52,53 could give 
an improved description. However, such functionals are 
still in an early stage of development and testing, and 
we prefer to focus first on exploring the extent to which 
conventional density functionals can reproduce experi- 
mental properties of systems in which orbital currents 
are present. 

Overall, significant progress has been made in the abil- 
ity to calculate the magnetoelectric coupling of real ma- 
terials in the context of density-functional theory. The 
methods described in Ref. 7 and 8 allow for the calcu- 
lation of both the electronic and lattice components of 
the spin (i.e., Zeeman) contribution to the magnetoelec- 
tric coupling. In principle at least, the lattice component 
of the orbital contribution could be computed using the 
methods of Ref. 54, while the remaining orbital electronic 
contributions can be computed from the formulas derived 
in Refs. 10 and 11 following the developments discussed 
here. We thus expect that the computation of all of the 
various contributions to the magnetoelectric coupling will 
soon be accessible to modern density-functional methods. 
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